APPENDIX A 



// 
// 
// 
// 
// 
// 
// 
// 
// 
// 
// 
// 



AUTHOR: A. CHRISTIAN TAHAN 
FIRST DRAFT: 02/10/00 



FILE: boolean. h 

FUNCTIONALITY: Boolean definitions and max min 

PROGRAM: required for all codes, checking for true/false for 



readability 



#ifndef _Boolean_type 

#define _Boolean_type 



#endif 

#ifndef FALSE 

#define FALSE 0 
#define TRUE 1 

#endif 

#endif 

#if !defined( MINMAX_DEFINED) 

#define MINMAX_DEFINED 

template <class T> T maxCT x, T y) 
{ 

return (x > y) ? x : y; 

}; 

template <class T> T min(T x, T y) 
{ 

return (x < y) ? x : y; 

}; 



#define AND && 
#define OR 1 1 
#define NOT ! 



#ifdef BOOL 



#define Boolean BOOL 



#else 



typedef int Boolean; 



#endif 



// 
// 
// 
// 

// I FILE: baseclas.cpp 

// I FUNCTIONALITY: test program for memory allocation (baseclas) 
// I PROGRAM: server space program 

// I COMMENTS: basic structure for database interaction including 
// I including file retrieval and bool program 

// I AUTHOR: A. CHRISTIAN TAHAN 
// I DATA FIRST VERSION: 02/10/00 
// 
// 
// 



#include "Baseclas. h" 
#include "Baseclas. hpp" 

class test 
{ 



public: 

t_real ma; 

static t_real num; 
#if (__BCPLUSPLUS__ == 0x340) 

char car[8]; 
#else 

char car[3]; 
#endif 

//this operator is always required if the class must be used in 
//ImpObjectList 

test & operator=(const test & a) 
{ 

ma=a . ma ; 

car[0]=a.car[0]; 
return C(*this)); 

} 

Boolean operator==(const test & a) 
{ 

return (ma==a.ma AND car[01==a.car[0]AND car[l]==a.car[l] 
AND car[2]==a.car[2] ); 

} 

testO 
{ 

carC0]=car[l]=car[2]='a' ; 



ma=num; 
num++; 

} 

test(int val) 
{ 

ma=val ; 
return; 

}; 

-testC) 
{ 

car[0]=car[l>car[2]='0' ; 
ma=0 . 0 ; 

} 

}; 

t_real test: :num=0; 



Boolean Test_Base_Class(t_index max_num) 
{ 

t_index i; 

ImpObjectList<test> mat; 
test compare; 

mwarn« n Base class Version " «BASECLAS_VERSION; 



mwarn«"testing base class allocating " «max_num«" elems of class 
test"; 

mat . Destroy _And_ReDimOnax_num) ; 

mwarn«"checking inizialization of " «max_num« M elems of class test"; 
for (i=0;i<max_num;i++) 
{ 

compare. ma=i+l; 
if (NOT (mat[i]==compare)) 
{ 

mwarn«"Inizialization error on element: "«i; 
return FALSE; 

} 

} 

mwarn«"checking access for " «max_num« , ' elems of class test"; 
for (i=0;i<max_num;i++) 
mat[i] .ma=i+l; 



for (i=0;i<max_num;i++) 



if (mat[i] .ma!=i+l) 
{ 

mwarn«"Inizialization error on element: 
return FALSE; 

} 

{ 

ImpSimpleTypeList<long> listl, Ust2 ; 
mwarn«"check operator ="; 
listl. Destroy_And_ReDim(70000UL); 
listl. Set(5); 
list2=listl; 

if (Ust2!=listl) 
{ 

mwarn«"operator = doesn't work on ImpSimpleList<long> 
return FALSE; 

} 

} 

mwarn«" check complete"; 
return TRUE; 

} 



// 
// 
// 
// 
// 
// 
// 
// 
// 
// 
// 
// 
// 



FILE: arraycla.cpp 
FUNCTIONALITY: array 
PROGRAM: required for all codes 

COMMENTS: to arrange the database without taking too much 
space 

AUTHOR: A. CHRISTIAN TAHAN 
DATE FIRST VERSION: 02/16/00 



#include "Arraycla.h" 
#include "Arraycla.hpp" 
#include "Diagclas.h" 
#include "Diagclas.hpp" 



MatrixOfDouble dummymat; 

Boolean Test_Array_Mat_Functions() 
{ 

MatrixOfDouble u,v,z; 
t_index dim=3,i,j; 

mwarn«"Array class Version: "«ARRAY_CLASS_VERSION; 

v . Destroy _And_ReDim(dim , dim) ; 
Check(v . DimC)==dim , 

"Error in the instruction \"v.Destroy_And_ReDim(dim,dim)\" 
v . Dim()="«v . DimO) ; 

u . Destroy _And_ReDim(dim , dim) ; 
u=v; 

for (i=0;i<dim;i++) 

for (j=0; j<dim; j++) 

Check(u[i][j]==v[i][j] /'Error in the instruction 

\"u=v\""); 

u.SetC2); 

for (i=0;i<dim;i++) 

for Cj=0; j<dim; j++) 

CheckCu[i][j]==2.0, 



ft 



"Error in the instruction \"u.Set(2)\" :u["«i«"] 

C"«j«"]="«u[i][j]); 



z=u*2.0; 

for (i=0;i<dim;i++) 

for (j=0; j<dim; 

Check(z[i][j>=4.0, 

"Error in the instruction \"z=u*2.0\" :z["«i«"] 

["«j«"]="«z[i][j]); 
u/=2.0; 

for (i=0;i<dim;i++) 

for (j=0;j<dim; 
Check(u[i][j]==1.0, 

"Error in the instruction \"u/=2.0\" :u["«i«"] 

="«u[i][j]); 

u*=2.0; 

for (i=0;i<dim;i++) 

for (j=0; j<dim; 
CheckCu[i][j]==2.0, 

"Error in the instruction \"u*=2.0\" :u["«i«"]["«j«"] 

="«u[i][j]); 

v=(u+2.0); 

for (i=0;i<dim;i++) 

for (j=0; j<dim; 
Check(v[i][j>=4.0, 

"Error in the instruction \"v=(u+2.0)\" :v["«i«"] 
C"«j«"]="«v[i][j]); 

v=v/2.0; 

for (i=0;i<dim;i++) 

for Cj=0; j<dim; 
Check(v[i][j]==2.0, 

"Error in the instruction \"v=v/2.0\" :v["«i«"]["«j«"] 

="«v[i][j]); 

u+=(v+z-u) ; 

for (i=0;i<dim;i++) 

for (j=0; j<dim; 
Check(u[i][j]==6.0, 

"Error in the instruction \"u+=(v+z-u)\" :u["«i«"] 
["«j«"]="«u[i][j]); 



u=u+v; 



for (i=0;i<dim;i++) 

for (j=0; j<dim; 
Check(u[i][j]==8.0, 

"Error in the instruction \"u=u+v\" :u["«i« ,, ]["«j«"] 

="«u[i][j]); 

u-=2.0; 

for (i=0;i<dim;i++) 

for (j=0; j<dim; 
Check(u[i][j]==6.0, 

"Error in the instruction \"u-=2.0\" :u["<<i«"]["«j«"] 

= n «u[i][j]); 

u-=v; 

for (i=0;i<dim;i++) 

for (j=0; j<dim; 
Check(u[i][j>=4.0, 

"Error in the instruction \"u-=v\" :u["«i«"]["«j«"] 

="«u[i][j]); 

u=ulv; 

for (i=0;i<dim;i++) 

for (j=0; j<dim; 

Check(u[i][;j>=24.0, 

"Error in the instruction \"u=ulv\" :u["«i«"]["«j«"] 

="«u[i][j]); 

VetDouble vet; 

vet . Destroy _And_ReDim(dim) ; 

vet.Set(2); 

vet=u I vet ; 

for (i=0;i<dim;i++) 

Check(vet[i]==144.0, 

"Error in the instruction \"vet=ul vet\" :vet["«i«"] 

="«vet[i]); 

Check (u!=v, "Error in boolean function"); 
u I =v ; 

for (i=0;i<dim;i++) 
f or( j=0 ; j<dim ; 

Check(u[i][j]=144.0, 

"Error in the instruction \"ul=v\" :u["«i«"] 
["«j«"]="«u[i][j]); 



MatrixOfDouble k, y, identity; 
t_real is_singular; 
t_index dim_k=4; 

'identity . Destroy _And_ReDim(dim_k , dim_k) ; 

identity[0][0]=1.0; 
identity[0][l]=0.0; 
identity[0][2]=0.0; 
identity[0][3]=0.0; 
identity[l][0]=0.0; 
identity[l] [1]=1.0; 
. identity[l][2]=0.0; 
identity[l][3]=0.0; 

identity[2][0]=0.0; 
identity[2] [1]=0.0; 
identity[2][2]=1.0; 
identity [2] [3>0.0; 
identity[3][0]=0.0; 
identity [3] [1]=0.0; 
identity[3][2>0.0; 
identity [3] [3] =1.0; 



k . Destroy_And_ReDim(dim_k , dim_k) ; 
y . Destroy _And_ReDim(dim_k , dim_k) ; 

k[0][0]=4.0; 

k[0][l>2.0; 
k[0][2]=2.4; 

k[0][3]=2.0; 

k[l][0]=0.0; 

k[l][l]=8.0; 

k[l][2]=3.6; 

k[l][3]=8.7; 

k[2][0]=5.1; 

k[2][l]=9.3; 

k[2][2]=2.9; 

k[2][3]=3.1; 

k[3][0]=7.23; 

k[3][l]=5.7; 

k[3][2]=1.9; 

k[3][3]=4.98; 



y=k; 



i s_si ngular=k . InverseO ; 

Check(is_singular!=0.0,' 'Routine InverseO doesn't work, 
is_singular="«is_singular); 

k=ylk; 

k.ChopO; 

Check(k==identity, "Routine InverseO doesn't work"); 

t_index dim_y=2; 
VetDouble vect; 

y . Destroy _And_ReDim(dim_y , dim_y) ; 

y[0][0]=2; 
y[0][i>i; 
y[i][0]=i; 
y[i][i]=2; 

y.EigenValues_And_EigenVectors(vect, k); 

MatrixOf Double eigval ,tr ,res; 

eigval . Destroy _And_ReDim(dim_y , dim_y) ; 

for (i=0;i<dim_y;i++) 
for (j=0; j<dim_y; 

if Ci-j) 
eigval[i][j]=vect[i]; 

else eigval[i][j]=0.0; 



tr ;Transpose_Of (k) ; 
res=tr I eigval Ik; 
res*=res; 

y*=y; 

k=res-y; 
k.ChopO; 

for (i=0;i<dim_y;i++) 
for (j=0; j<dim_y; 

Check(k[i][j]<=PRECISION, "routines 
EigenValues_And_EigenVectors don't work"); 



DiagMatrixOfDouble diag; 
Mat rixOf Double kk; 
t_index dim_diag=4; 

diag. Destroy _And_ReDim(dim_diag, dim_diag) ; 

diag[0][0]=4.0; 
diag[l][l>6.0; 
diag[2][2]=5.1; 
diag[33 [33=7.23; 

k«=diag; 

y . Destroy_And_ReDim(dim_diag , dim_diag) ; 

y[0][0]=4.0; 

y[0] [13=0.0; 
y[0][2]=0.0; 

y[0][3]=0.0; 

y[l] [0]=0.0; 

y[13[13=6.0; 

y[i3[23=0.0; 

y [13 [33=0.0; 

y [23 [03=0.0; 

y[23[i3=0.0; 

y [23 [23=5.1; 

y[23[33=0.0; 
y [33 [03=0.0; 
y[33[i3=0.0; 
y[33[23=0.0; 
y[33[33=7.23; 

Check(y==k, "Routine Charige_Diag2Full don't work"); 
return TRUE; 

} 



// Housholder method transforms a symmetric matrix into a tridiagonal one 
void Mat rixOf Double: :Householder() 
{ 

VetDouble x, y, u; 
MatrixOfDouble H, U, UT, I; 



t_index dim, i, k; 



t_real sum; 



dim=(*this)[0].DimO; 

// B is a symmetric matrix of order greater than two 
// the symmetric property is not tested 
Assert(dim>2); 

x . Destroy _And_ReDim(dim) ; 
y . Dest roy_And_iReDim(dim) ; 

H . Dest roy_And_ReDim(dim , dim) ; 
U . Dest roy_And_ReDim(dim , dim) ; 
UT . Dest roy_And_ReDim(dim , dim) ; 

I . Dest roy_And_ReDim(dim , dim) ; 

f or(i=0 ; i<dim ; i++) 
I[i][i>1.0; 

f or(k=0 ; k<dim-2 ; k++) 
{ 

UT.Set(0); 

H.Set(0); 

y.Set(0); 

for(i=0;i<dim;i++) 

x[i]<*this)[i][k]; 

for(i=0;i<=k;i++) 
y[i]=x[i]; 

sum=0; 

for(i=k+l;i<dim;i++) 

sum+=pow(x[i] ,2); 
sum=sqrt(sum) ; 
if(x[k+l]>0) 

sum=-sum; 

y[k+l]=sum; 
UT[0]=x-y; 

// calculate the x-y norm 
sum=0; 

for(i=0 ; i<dim ; i++) 

sum+-pow(UT[0] [i] ,2) ; 
sum=sqrt(sum); 

UT[0]/=sum; 



U.Transpose_Of(UT); 



U=UIUT*2; 
H=I-U; 

(*this)=HI(*this)IH; 
} 

return; 
} 



#define SIGN(a,b) C(b)<0 ? -fabs(a) : fabs(a)) 

void MatrixOf Double: : Eigenvalues J\nd_EigenVectors(VetDouble& eigenvalues, 

Mat rixOf Double& 

eigenvectors) const 
{ 

t_signed i, I; 

t_index m, j, k, iter, dim; 

t_real s, r, p, g, f, dd, c, b; 

VetDouble external, diagonal, eig_values; 

MatrixOfDouble tridiagonal; 

tridiagonal = (*this); 

dim= (*this).Dim(); 

eigenvalues . Destroy _And_ReDim(dim) ; 

external . Destroy _And_ReDim(dim) ; 

diagonal. Destroy _And_ReDim(dim); 

ei genvectors . Destroy _And_ReDim(dim , dim) ; 

ei genvectors . Set(0 . 0) ; 

tridiagonal , HouseholderO ; 

f or(i=l ; i<(t_signed)dim ; i++) 
{ 

external [i]=tridiagonal [i -1] [i] ; 
diagonal [i]=tridiagonal [i] [i] ; 
eigenvectors[i][i]=1.0; 
} 

f or(i=l ; i <(t_signed)dim ; i++) 

external [i -l]=external [i] ; 

external [dim-l]=0.0; 



for (l=0;l<(t_signed)dim;l++) 
{ 

iter=0; 
do { 

for (m=l;m<dim-l;m++) 
{ 

dd=fabs(diagonal [m])+f abs(diagonal [m+1] ) ; 
if (fabs(external[m])+dd == dd) 
break; 

} 

if ((t_signed)m != I) 
{ 

if (iter++ == 30) 

cout«"Too many iterations"; 
g=(diagonal [1+1] -diagonal [I] )/(2 . 0*external 



Ci]); 

+SIGN(r,g)); 



r=sqrt(Cg*g)+l.0); 

g=diagonal [m] -diagonal [l]+external [I] /(g 

s=c=1.0; 
p=0.0; 

for (i=m-l;i>=l;i--) 
{ 

f=s*external[i] ; 
b=c*external[i] ; 
if (fabs(f) >= fabs(g)) 
{ 

c-g/f; 

r=sqrt((c*c)+l . 0) ; 
external [i+l]=f *r ; 
c *= (s=1.0/r); 
} 

else 

{ 

s=f/g; 

r=sqrt((s*s)+1.0); 
external [i+l]=g*r ; 
s *= (c=1.0/r); 
} 

g=diagonal [i+1] -p ; 
r=(diagonal [i] -g)*s+2 . 0*c*b ; 
p=s*r; 

diagonal [i+l]=g+p ; 
g=c*r-b; 

for Ck=0;k<dim;k++) 



{ 

f=eigenvectors[k] [i+1] ; 
eigenvectors[k] [i+1] 

=s*eigenvectors[k][i]+c*f ; 

eigenvectors[k] [i] 

=c*eigenvectors[k] [i]-s*f ; 

} 

} 

diagonal [l]=diagonal [1] -p ; 
external [l]=g; 
external [m] =0.0; 
} 

} while ((t_signed)m != 1); 

} 

// transposition of eigenvectors matrix in order to have 
autovectors on the rows 

// and not on the columns 
t_real temp; 

f or(i=0 ; i<(t_signed)dim-l ; i++) 
f or( j=i+l ; j<=dim-l ; 
{ 

temp = eigenvectors[i][j] ; 
eigenvectors[i][j] = eigenvectors[j][i] ; 
eigenvectors[j][i] = temp; 
} 

return; 

}; 



yy ******************************* 

// * 

// Matrix of double 

// * 

****************************************************************** 



// Perform a low up decomposition from a square matrix 
void MatrixOf Double: :Low_Up_Dcmp(MatrixOf Double &L, MatrixOfDouble &U) 
{ 

t_index i,j,k,n; 
t_real sum; 

Assert((*thi s) . Dim_Row(>=(*this) . Dim.ColQ) ; 



n=(*this).Dim_RowC); 



L . Destroy_And_ReDim(n ,n) ; 
U . Destroy_And_ReDim(n ,n) ; 

for(j=0;j<n; 

L[j][j]=1.0; 

j=0; 

for(i=0;i<n;i++) 

U[j][i]=(*this)[j][i]; 

for(i=l;i<n;i++) 

L[i][j]=C*this)[i][j]/U[j][j]; 

for(j=l;j<n;j++) 
{ 

for(i=j ;i<n;i++) 
{ 

sum=0 . 0 ; 

for(k=0;k<=j-l;k++) 

sum+=L[j][k]*U[k][i]; 
U[j][i]=(*this)[j][i]-sum; 
} 

for(i=j+l;i<n;i++) 
{ 

sum=0.0; 

for(k=0;k<=j-l;k++) 

sum+=L[i][k]*U[k][j]; 
L[i][j]=((*this)Ci][j]-sum)/U[j][j]; 
} 

} 

return ; 
} 

// Solve a linear equation Ax=b where A is a lower triangular or 
// upper triangular matrix 

void MatrixOf Double: :Solve_Triangular(VetDouble &y, VetDouble &b) 
{ 

t_index i,k,n; 
t_real sum; 

n=b.Dim(); 

//Assert(tri_mat . Dim_Row()==tri_mat . Dim_ColO) ; 
Assert(n==(*this) . Dim_Row()) ; 



y . Destroy_And_ReDim(n) ; 

if((*this)[0][l>=0.0) 
{ 

y[0]=b[0]; 
for(i=l;i<n;i++j 
{ 

sum=0.0; 

for(k=0;k<=i-l;k++) 

sum+=(*this)[i][k]*y[k] ; 

y[i]=b[i]-sum; 
} 

} 

else{ 

y[n-l]=b[n-l]/(*this)[n-l][n-l]; 
for(i=n-l;i>0;i--) 
{ 

sum=0.0; 

for(k=n-l;k>=i ;k--) 

sum+=(*this)[i-l] [k]*y[k] ; 

y[i-l]=(b[i-l]-sum)/(*this)[i-l][i-l]; 
} 

} 

return; 
} 

void Matri xOf Double : :Transpose_Of (const MatrixOfDouble & x) 
{ 

Destroy_And_ReDim(x. Dim_Col() ,x . Dim_Row()) ; 
t_index i,j; 

f or(i=0 ; i<x . Dim_Row() ; 

for( j=0 ; j<x . Dim_ColO ; 

C*this)[j][i]=x[i][j]; 

return; 
} 

t_real MatrixOfDouble: :Inverse() 
{ 

MatrixOfDouble L,U,B,X; 
VetDouble y; 
t_index N,k; 



t_real det; 

N=(*this).Dim_ColO; 
Assert((*this).Dim_Row()==N); 

(*this).Low_Up_Dcmp(L, U); 

B . Destroy _And_ReDim(N , N) ; 
X . Destroy _And_ReDim(N , N) ; 

det=1.0; 

for(k=0;k<N;k++) 
{ 

B DODO =1.0; 

det*=U[k][k]; 

} 

if(fabs(det)<=10E-10) 

mwarn«"Singular matrix in low up decomposition"; 

for(k=0;k<N;k++) 
{ 

L.Solve_Triangular(y, B[k]); 
U . Solve_Triangular(X[k] , y) ; 
} 

(*this) .Transpose_Of (X) ; 

return det; 
} 



// 
// 

// 

// I 

// I 

// I FILE: diagclas.cpp 

// I FUNCTIONALITY: diagonal matrix 

// I PROGRAM: required to all codes 

// I COMMENTS: header template classes, for access to header files 

// I AUTHOR: A. CHRISTIAN TAHAN 

// I DATE FIRST VERSION: 02/13/00 

// | 

// _ 



#include "Diagnost.h" 
#include "Diagclas.hpp" 
#include "Arraycla.h" 
#include "Arraycla.hpp" 

Boolean Test_Math_Diagonal_Matrix_Function O 
{ 

DiagMatrixOfDouble u,v,z,c; 
t_index dim=10, i ; 

mwarn«"diagclas version"«DIAG_CLASS_VERSION; 

u . Destroy _And_ReDim(dim , dim) ; 
v . Destroy _And_ReDim(dim , dim) ; 
z . Dest roy_And_ReDim(dim , dim) ; 
c . Destroy _And_ReDim(dim , dim) ; 

Check (v.Dim0==dim, "Error in Destroy_And_ReDim() function"); 

u.Set(Z); 
v.Set(3); 
z.Set(2); 
z*=u*v*u*u*5; 

c.Set(240); 

Check(c==z, "diagonal matrix routines for the product don't work"); 

c+=u+80+v ; 
z.Set(325); 



Check(c==z, "diagonal matrix routines for the addition don't work"); 
Check(z[dim-l][dim-l]==c[0][0], ,, diagonal matrix routine operator [] don't 
work"); 

c/=25; 

for (i=0;i<dim;i++) 
Check(c[i][i]==13.0, "diagonal routines operator/= don't work"); 

c=ul vl v; 

for (i=0;i<dim;i++) 

Check(c[i][i]==18.0, "diagonal routines operator I don't work"); 

c=ul v; 

for (i=0;i<dim;i++) 
Check(c[i][i]==6.0, "diagonal routines operator! don't work"); 

VetDouble vect, ris; 

vect . Destroy _And_ReDim(dim) ; 
vect.Set(2); 

ris=cl vect; 

for (i=0;i<dim;i++) 

Check(ris[i]==12.0, "diagonal routines operator I don't work"); 

v. Set (3); 

cl=v; 

for (i=0;i<dim;i++) 

Check(c[i][i]==18.0, "diagonal routines operatorl= don't work"); 
u[0] [0]=12.0; 
u[l][l]=22.5; 
u[2][2]=343.1; 
u[3][3]=125.4; 
u[4][4]=2.0; 
u[5] [5] =458.1; 
u[6] [6] =75. 2; 
u[7][7>45.126; 
u[8] [8] =75. 2; 
u[9] [9] =45. 3; 

v=u ; 

u.Inverse(); 
c=v I u ; 
c.Chop(); 



f or(i=0 ; i<dim ; i++) 
Check(c[i][i>=1.0, 

"diagonal matrix routines InverseO don't work : c["«i«"]="«c 

Ci][i]); 
DiagMatrixOf Double eigvect; 
v.Set(3); 

v . Ei genValues_And_EigenVectors(vect , eigvect) ; 

for(i=0;i<dim;i++) 
u[i][i]=vect[i]; 

Check(u==v,"error in function EigenValues_And_EigenVectors H ); 

return TRUE; 

} 



